A&A manuscript no. 

(will be inserted by hand later] 



Your thesaurus codes are: 
02(12.07.1; 03.13.1; 03.13.4) 



ASTRONOMY 

AND 
ASTROPHYSICS 



A fast direct method of mass reconstruction for gravitational 
lenses 






M. Lombardi^'^ and G. Bertin^ 



^ European Southern Observatory, Karl-Schwarzschild Strafie 2, D 85748 Garching bei Miinchen, Germany 
^ Scuola Normale Superiore, Piazza dei Cavalieri 7, I 56126 Pisa, Italy 



Received 25 March 1999; accepted 17 May 1999 



1—5 



> 



o 



Oh' 

I ' 

o- 



X 



Abstract. Statistical analyses of observed galaxy distor- 
tions are often used to reconstruct the mass distribu- 
tion of an intervening cluster responsible for gravitational 
lensing. In current projects, distortions of thousands of 
source galaxies have to be handled efhciently; much larger 
data bases and more massive investigations are envisaged 
for new major observational initiatives. In this article we 
present an efficient mass reconstruction procedure, a di- 
rect method that solves a variational principle noted in an 
earlier paper, which, for rectangular fields, turns out to re- 
duce the relevant execution time by a factor from 100 to 
1000 with respect to the fastest methods currently used, 
so that for grid numbers N = 400 the required CPU time 
on a good workstation can be kept within the order of 1 
second. The acquired speed also opens the way to some 
long-term projects based on simulated observations (ad- 
dressing statistical or cosmological questions) that would 
be, at present, practically not viable for intrinsically slow 
reconstruction methods. 
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1. Introduction 

In the context of weak or statistical lensing, the problem of 
the determination of the dimensionless mass density distri- 
bution k(6) from a map of the reduced shear g{9) has been 
considered in detail by various authors, using either simu- 
lations or analytical calculations (e.g., Bartelmann 1995, 
Schneider 1995, Seitz & Schneider 1996, Squires & Kaiser 
1996, Lombardi & Berlin 1998a, 1998b). 

The mass inversion is usually performed starting from 
the vector field u{6) defined in terms of the measured re- 
duced shear g{0) (Kaiser 1995). In the ideal case where 
the measured shear g(6) is just the true shear go{9), the 
vector field uo{9) can be shown to satisfy the relation 



«o(e)=Vln[l-Ko(0)] =Vko(0), 
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(1) 



where kq is the true dimensionless mass map and ko{9) = 
ln[l — K,o{9)\. However, because of statistical and mea- 
surement errors, u is not necessarily curl-free, and thus 
K can be determined only approximately. In a separate 
paper (Lombardi & Berlin 1998b) we have shown that 

— The statistical errors on k{9) are minimized if this 
function is calculated as 



k{9) = K+ / H^^{9,9')-u{9')d^e' 



(2) 



where k is a constant (introduced to take into account 
the sheet invariance), H^^{9,9') is the noise-filtering 
kernel (Seitz & Schneider 1996), and H, is the field of 
observation. 
— The same mass map can be obtained by solving the 
equations 



Vk ■ n ~ u ■ n 



on 9il 



(3) 
(4) 



where n is the unit vector perpendicular to the bound- 
ary dVt of the field of observation $1. Hence, the kernel 
H^^{9,9') can be identified as the Green function of 
this Neumann boundary problem. 
Equations (||) and (||) are precisely the Euler equations 
associated with the functional 



5* = 



1 



\Vfii9)-u{9)\\ d^e 



(5) 



In other words, the functional S is minimized when 
k{9) is calculated using Eq. (g) or, equivalently, by 
solving Eqs. (||) and (||). 

To these three formulations of the mass inversion problem 
there correspond three practical methods to calculate k{9) 
from a given set of data. 

— The first method considered is based on a direct cal- 
culation of the kernel H^^ (Seitz & Schneider 1996). 
Once this kernel has been calculated for a given field 
Q, the mass inversion is straightforward (note that 
the kernel H^^ depends on the field of observation). 
A problem with this method is that a calculation of 
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H is expensive in terms of memory requirements and 
computation time. In fact, in order to compute /t on a 
square grid of A^ x A^ points, H^^ must be calculated 
on a multidimensional grid ofNxNxNxN points. 
Moreover, 2N^ multiplications are needed to evaluate 
Eq. (|), and thus the method is of order 0{N'^). Be- 
cause of the large memory needed to allocate H^^, cal- 
culations can be performed only with a limited value 
of A^ (typically A^- 50). 

— The introduction of a method that directly solves the 
Neumann problem allows one to go beyond many of the 
limitations of the H^^ method. Equations (||) and (^ 
can be solved using an over-relaxation method (Seitz 
& Schneider 1998). In this case k{9) is calculated di- 
rectly, and thus we need to allocate only N x N real 
numbers. Moreover, the method can be applied with- 
out difficulties to "strange" geometries il (while the 
previous method is straightforward only when applied 
to rectangular or circular fields). The over-relaxation 
method is quicker than the kernel method, being ap- 
proximately of order 0[N^). 

— A direct method to minimize the functional (0) will be 
presented in this paper. As we will see, this method 
has several advantages and turns out to be very ef- 
ficient from a computational point of view, being of 
order 0[N^\ogNy Moreover, it is extremely easy to 
implement (two implementations for rectangular fields 
fl written in C and in IDL are freely available on re- 
quest). 

We should stress that, as proved in an earlier paper 
(Lombardi & Bertin 1998b), the three formulations are 
mathematically equivalent. Thus it would not be surpris- 
ing to find that proper numerical implementations per- 
form, for large values of the grid number A^, in a similar 
manner as far as accuracy and reliability arc concerned. 
In practice, for finite values of A^, the third method turns 
out to be characterized by small errors, often smaller than 
those associated with the other two procedures. 

2. A direct method to solve the variational 
principle 

Direct methods in variational problems are well-known 
especially in applied mathematics (see, e.g., Gelfand & 
Fomin 1963). Suppose that one can find a complete set 
of functions {fa} on the domain fi (the full definition of 
"complete" will be given below), so that any function on 
fl can be represented as a linear combination of the form 



k{e) 



a=l 



.Ue) 



(6) 



More precisely, we assume that for any function ii{6), 
there is a choice for the coefficients {ca} such that 



/ 



Vk(6») 



a=l 



caVUie) 



d^0 = O 



(7) 



Let us now introduce a sequence of trial mass maps 

71 



Q = l 



We further require that the function kI"1 minimizes the 
functional S: in other words, Ci ,C2 , ■ ■ ■ ,Cn are chosen 
so that the functional S has minimum value. This obvi- 
ously happens when 



dS 







for a = 1,2,. 



(9) 



Solving this set of n equations, we obtain the n coeffi- 
cients Cq , and thus the function k^"!. By repeating this 
operation for a sequence of values of n, we find a sequence 
of functions k^"'. These functions, under suitable assump- 
tions (verified in our problem) , have the following proper- 
ties (see Gelfand & Fomin 1963 for a detailed discussion): 
(i) Let us call S'^"' the value of S when k is replaced by 
the function «;["! . Then, obviously, the sequence S*'"! is not 
increasing, (ii) If the set {fa} is complete, then the func- 
tions ^["1 converge to the solution k of the problem. This 
method thus provides a way to obtain the function k{9) 
with desired accuracy. 

The method described here can be easily applied to 
our problem. In fact, by expanding k^^\d) as in Eq. (||), 
we find 



dS 



= /v/„(0) 

Jn 



J2c^;^vf,{e)-u{e) 

73=1 



d26' = 0.(10) 



The previous equation, for a — 1,2, ... ,n, represents a 
linear system of n equations for the n variables {cq }. Its 
solution is thus the set of coefficients to be used in Eq. (||). 
However, we note that care must be taken in the choice 
of a complete set of functions. Let us define, for the pur- 
pose, the product {v, w) between two generic vector fields 
v{6) and w{9) as 



Jn 



(11) 



As our problem involves Vk, the completeness has to be 
referred to the set of the gradients. In other words, the set 
{fa} is complete if 



{Vfa,Vk)=0 



(12) 



for every a implies k{9) — constant. It is easy to show 
that this condition is equivalent to Eq. (0). 

The direct method can be further simplified if a set of 
functions {fa} can be taken to satisfy a suitable orthonor- 
mality condition, so that the gradients of the functions 
{fa} satisfy 



{^fa,^fp)^5ap 



(13) 
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where Sap = 1 for a = /? and otherwise. Then Eq. (|To|) 
can be rewritten simply as 



(V/„,«> 



(14) 



Thus, with the use of an orthonormal set of functions 
we have secured two important advantages: (i) The lin- 
ear system (noh has been diagonalized, so that its solution 
is trivial, (ii) The coefficients Cq no longer depend on n: 
that is, the coefficients of the exact solution are given by 
Ca = (V/q,m). 

Because of these advantages, whenever possible an or- 
thonormal set of functions should be used. We note, how- 
ever, that the orthonormality condition (|l3| ) depends on 
the field of observation fi. Even if the existence of an or- 
thonormal set of functions is always guaranteed by the 
spectral theory for the Laplace operator (see Brezis 1987), 
for "strange" geometries, it may be non trivial to find a 
complete orthonormal set of functions. In such cases, we 
need to solve the linear system ([lO|). 

The direct method described above has several ad- 
vantages with respect to the "kernel" method and to the 
over-relaxation method: (i) The method is fast in the case 
where an orthonormal set of functions can be found. In 
fact, we need only to evaluate one integral for each coeffi- 
cient Ca that we want to calculate, (ii) The method does 
not require a large amount of memory: we need to retain 
only the n values of the coefficients Ca- (iii) The precision 
of the inversion is driven in a natural way by the value of 
n. Typically, the larger a is, the smaller the length scale 
of fa (see below), (iv) In some cases, the decomposition 
of the mass density k(6) in terms of the functions /„ can 
be useful. 

3. Rectangular fields 

When the field fl is rectangular, an orthonormal set of 
functions can be written easily. Here we consider the spe- 
cial case when fi is a square of length ir (in some suitable 
units); any rectangular field can be handled in a simi- 
lar manner. In the case considered, an orthonormal set of 
functions is given by 



/q/3 (d) — riafi cos aOi cos (592 



(15) 



with (a, /3) e N'^ \ (0, 0). The normalization Uap is defined 
as 

^ for a = or /3 = , 

(16) 



.., = <! '^v/^ 



7rVa2+752 



otherwise . 



The function /oo is not defined. Note that here we use two 
indices for the set. Cosines must be used in order to have 
a complete set (see Eqs. (0), (|l^), and Appendix A). Our 
problem is solved in terms of the coefficients Cap'- 



^ap 



T-ap I [aui(0) sina^i COS/ 



f3u2 (d) COS aOi sin (362] d^ 
'^W — /^Cc(/3riQ^ cosa6'i cos/36'2 . 



(17) 
(18) 



ap 



We now observe that the particular choice of the orthonor- 
mal set {/a/3} allows us to use fast Fourier transform 
(FFT) techniques to evaluate Eqs. (^ and (^. The use 
of FFT makes the direct method very efficient: in particu- 
lar the method becomes of order 0(iV^ log N) . Moreover, 
several optimized FFT libraries are available. 

The optimal truncation for the series dl8) is deter- 
mined by the adopted grid numbers: for a grid oi N x M 
points, a should run from to iV — 1, and f3 from to 
M — 1 (this is standard practice for FFT libraries). 

4. Performance 

Our method has been implemented in C and in IDL. The C 
version uses the library FFTW ( "Fastest Fourier Transform 
in the West," version 2.0.1) to perform discrete Fourier 
transforms (DFT). This library, written by Matteo Frigo 
and Steven G. Johnson, is considered the quickest DFT 
library publicly available. The performance of our di- 
rect method is compared with that of the over-relaxation 
method, also implemented in C. The procedure used in the 
tests is summarized in the following points: 

1. A simple model for the dimensionless mass distribu- 
tion has been chosen. Then the mass distribution kq is 
calculated on a grid of iV x TV points. 

2. The associated field Uq is calculated on the same grid 
using a 3-point Lagrangian interpolation in order to 
numerically evaluate the derivatives that are needed. 

3. Noise is added to the vector field Uq using an analytical 
model for the noise derived earlier (Lombardi & Bertin 
1998b). In practice, the various Fourier components 
of the noise are added using a suitable model for the 
power spectrum. 

4. The resulting noisy u map is inverted using the over- 
relaxation method and the present direct method. The 
two dimensionless mass maps obtained are then com- 
pared. Moreover, the inversion times are recorded. 

The results obtained in the tests are the following: 

— The two mass densities obtained are consistent with 
each other. 

— Because of the set of functions used, the errors pro- 
duced by the direct method are larger on the bound- 
ary of the field. For this reason, we suggest that a one 
pixel strip around the field should be discarded. The 
area discarded is very small. 

— Some tests have been performed by providing Uq to the 
inversion procedures. This allows us to compare the 
reconstructed mass density with the original map kq. 
From such tests we have noted that the discretization 
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Fig. 1. Execution time per call vs. grid number N. The 
solid line refers to the direct method applied to "good 
numbers" (values 2iV— 1 that can be factorized with small 
primes), the long-dashed refers to the direct method ap- 
plied to "bad numbers" (27V — 1 prime), and the short- 
dashed line to the over-relaxation method. 



errors of the direct method are slightly smaller than the 
discretization errors of the over-relaxation method. 
~ The results of the two methods differ because of the 
sheet invariance: in particular, the direct method al- 
ways gives a "reduced" mass map k with vanishing 
total mass. 

Regarding the second item, we note that the error is re- 
lated to the finite sampling scale of the method; the er- 
ror affects only the outermost pixel because of the proper 
choice of the truncation (see comment at the end of 
Sect. 3). 

The measured execution times are plotted in Fig. 1 for 
different values of TV. These are the averaged CPU execu- 
tion times for a single reconstruction on a SUN Ultra 1 
workstation. From this figure it is clear that the direct 
method is much faster than the over-relaxation method. 
Here we should recall that, because of some character- 
istics of the FFTW library, the execution time of the di- 
rect method can change significantly even for neighbour- 
ing values of N. In particular, the inversion is faster when 
(2N — 1) can be factorized with small prime numbers, and 
is slower in other cases (see Fig. 1). For example, the ex- 
ecution time (on a SUN Ultra 1) changes from 2.942 to 
0.232 seconds when N changes from 121 to 122. Finally, 
we observe that our implementation of the direct method 
is not optimal: in fact, with a different (non-trivial) use 
of FFT one might gain an additional factor of 4 on the 
execution time. 

Besides the appealing aspects of simplicity inherent to 
the direct method described in this paper, we should note 
that gaining three orders of magnitude in CPU time will 
make it possible to undertake a few long-term projects 
of simulated observations (in particular, with the goal of 
a statistically sound investigation of the quality of mass 




Fig. 2. A typical result of mass reconstruction; at the 
adopted distance for the lensing cluster, the side of the 
square field, 10 arcminutes, corresponds to approximately 
2.88 Mpc. From top to bottom, true dimensionless mass 
distribution, reconstructed distribution (from the direct 
method), and difference between maps of the variable k 
derived from direct and over-relaxation methods. The very 
small residuals show that the two methods are practically 
equivalent in terms of accuracy. 

reconstruction; but other objectives might be formulated, 
e.g. in the cosmological context) that would remain prac- 
tically out of reach for other intrinsically slow reconstruc- 
tion methods. 

5. Examples of simulated reconstructions 

In addition to the reconstructions from "synthetic" data 
as described in the previous Section, we have performed 
several additional tests in order to demonstrate the relia- 
bility of our method. The tests, designed with the aim 
to reproduce the main features of a "real" reconstruc- 
tion, have been made following a straightforward proce- 
dure (see, e.g., Lombardi & Bertin 1999 for similar simu- 
lations). 
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First we have generated a population of source galaxies 
using a pseudo-random number generator. Here a source 
galaxy is represented by its position and by its ellipticity 
(see Eq. (11) of Lombardi & Bertin 1999, following Seitz & 
Schneider 1997). Positions are drawn from a homogeneous 
distribution (with a density of 70 galaxies per square ar- 
cmin), while ellipticities are drawn from a truncated Gaus- 
sian distribution with variance a"^ = (0.3)^. Sources are 
assumed to have all the same redshift Zs — 1.5. Source el- 
lipticities are then transformed into observed ellipticities. 
For simplicity, the observed galaxy positions are assumed 
to be equal to the source positions: in other words, no 
depletion effects are included in the simulations. 

Then the calculation of the observed ellipticities has 
been done by referring to a cluster of galaxies placed 
at Zd = 0.3 with total mass inside the 10' x 10' field 
1.5 X 10^^ M0. For the purpose of introducing the lens- 
ing effects, we only need to specify the dimensionless pro- 
jected mass map k,{0). For simplicity, we have used a den- 
sity distribution made of three symmetrical components; 
each component is described by the analytical model out- 
lined by Schneider et al. (1992, p. 244), which, at large 
radii, is approximately isothermal. 

By averaging the observed ellipticities, we have then 
obtained a map of the reduced shear g{0) and, from that 
map, the vector field u{6). The mass inversion has been 
performed using the direct method and the over-relaxation 
method. The two mass distributions have been then com- 
pared. 

An example of typical results is shown in Fig. 2. Here, 
from top to bottom, we display the original cluster mass 
distribution, the reconstruction obtained using the direct 
method, and the residuals, i.e. the difference between the 
reconstructed maps from the direct method and from the 
over-relaxation method. As the figure clearly shows, dif- 
ferences are mainly confined to the boundary of the field 
where they are found to be of the order of 0.0002, well 
below the statistical errors of the reconstruction. In the 
inner field the differences are about one order of magni- 
tude smaller. Note also that the wavy overall appearance 
of the reconstructed map is normal for weak lensing re- 
constructions, resulting from the relatively low number of 
source galaxies involved (see Lombardi & Bertin 1998b for 
a discussion of the statistical aspects of the problem). 
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Appendix A: Completeness of {fap} 

In this appendix we will verify explicitly that the set of 
functions defined in Eq. ( [l5[ ) is complete, in the sense that 
Eqs. (y) and (Q) can be recovered. For the purpose, we 



will apply the Fourier theorem (Brezis 1987). In the fol- 
lowing, we will assume that u{0) is a smooth vector field 
(we stress that this condition is needed only for the proof 
provided below; the method remains applicable to more 
general cases). 

Let us consider a solution of the form (@). Then, 
because of the orthonormality condition dl3), we have 
(V/q, Vk) =Ca = (V/q, it), so that 

(V/a, Vk -u) = I V/a • (Vk - m) d^e* = . (A.l) 

Now we observe that the previous equation holds for 
any a: then it holds also for any linear combination 

/^Ea^-/" of {/a}. Thus 

== I vf -{vk- u)A^e = - / /v • (Vk - m) d^e 



f{Wk- 



nde 



(A.2) 



90 



In the last step we have integrated by parts (n is the unit 
vector orthogonal to the boundary Qfi of fi). 

We now use this equation to show that the chosen set 
of functions, described by Eq. (p^, is complete, while, e.g., 
a similar set made of sine functions would not be complete. 
By the nature of the chosen set of functions we already 
know that we can properly represent any smooth function 
/. Using this property, we want to show that the two terms 
V • (Vk — u) and (Vk — u) ■ n entering in the r.h.s. of 
Eq. (A.2) vanish on fi and on d^ respectively. 

For the purpose, we observe that if cosines are used 
as set of functions, we can "build" any function / pro- 
vided that the function has periodic derivatives on the 
boundary. In particular, if A C 51 is an arbitrary open 
subset of fi, there is a function / that is positive on A and 
vanishes on $7 \ ^4. Now suppose per absurdum that the 
solution obtained from the direct method does not satisfy 
Eq. (H), so that, e.g., V • (Vk - it) > on a point 6* £ fi. 
Then, for the sign persistence theorem, this quantity must 
be strictly positive in a neighborhood A oi 9* . However, 
if we take a function / which is positive on A and van- 
ishes elsewhere, the rhs of Eq. ( |A.2|) will be positive, while 
the Ihs vanishes, which is contradictory. This proves that 
Eq. (H) is verified by cosines. 

In a similar manner, now that we have "disposed of" 
the first term, we observe that using cosines we can build 
a function / that vanishes everywhere on the boundary 
of dfl except for a neighborhood. In other words, given 
an open subset B C dil of the boundary 957, there is a 
function / that is positive on B and vanishes on dfl \ B. 
Note that this property would not be satisfied if the set of 
functions {/„} were based on sines. Using a proof similar 
to the one given above, we obtain that (Vk — it) • n must 
vanish, thus leading to Eq. (0). 

One might worry that, on the boundary, the chosen 
set of functions of Eq. ( [l5|) has zero derivative in the di- 
rection of n, i.e. V/q/3 • n = 0. This might suggest that 
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the boundary condition (m cannot be reproduced. In re- 
ahty, although this is true pointwise, this does not affect 
the convergence in C^ that is relevant for our problem (see 
Eqs. (M) and (112)). This important point is best illustrated 
by the following example. 

Suppose that we measure a constant field u{0) = (1, 0) 
on a square field Q of side n. The obvious solution for k 
in this case is k{9) = 9i. Suppose now that we try to 
use a set of functions made of sines. Then the correspond- 
ing coefficients Caf3 would be proportional to the integrals 
a J^ 01 cos a9i d9i J^ sin/36'2 d02 — 0. Hence, every coeffi- 
cient Cap vanishes. This proves that a set based on sines 
is not complete (condition (|2|) is not satisfied or, equiva- 
lently, a curl-free vector field u leads to a vanishing mass 
map). On the other hand, if we use the set ([l5|), the co- 
efficients Cap do not vanish and the corresponding mass 
distribution is given by 



k(e) = - V — ^cosa6'i. (A.3) 

a odd 

This function can be shown to reduce to k{6) = 9i — tt/2. 
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